/* Integration points (quadratures)
   Jukka Toivanen, Jan Stebel
*/

#include "quadpoint.h"




template<>
void QuadPoint<3>::setGSF()
{
  GSF(0,0) = -3 + 4*(ksi[0]+ksi[1]+ksi[2]);
  GSF(1,0) = -3 + 4*(ksi[0]+ksi[1]+ksi[2]);
  GSF(2,0) = -3 + 4*(ksi[0]+ksi[1]+ksi[2]);
  GSF(0,1) = -1 + 4*ksi[0];
  GSF(1,1) = 0;
  GSF(2,1) = 0;
  GSF(0,2) = 0;
  GSF(1,2) = -1 + 4*ksi[1];
  GSF(2,2) = 0;
  GSF(0,3) = 0;
  GSF(1,3) = 0;
  GSF(2,3) = -1 + 4*ksi[2];
  GSF(0,4) = 4*(1-2*ksi[0]-ksi[1]-ksi[2]);
  GSF(1,4) = -4*ksi[0];
  GSF(2,4) = -4*ksi[0];
  GSF(0,5) = 4*ksi[1];
  GSF(1,5) = 4*ksi[0];
  GSF(2,5) = 0;
  GSF(0,6) = -4*ksi[1];
  GSF(1,6) = 4*(1-ksi[0]-2*ksi[1]-ksi[2]);
  GSF(2,6) = -4*ksi[1];
  GSF(0,7) = -4*ksi[2];
  GSF(1,7) = -4*ksi[2];
  GSF(2,7) = 4*(1-ksi[0]-ksi[1]-2*ksi[2]);
  GSF(0,8) = 0;
  GSF(1,8) = 4*ksi[2];
  GSF(2,8) = 4*ksi[1];
  GSF(0,9) = 4*ksi[2];
  GSF(1,9) = 0;
  GSF(2,9) = 4*ksi[0];
}


template<>
void QuadPoint<3>::setSF()
{
  SF(0,0) = 1-3*ksi[0]-3*ksi[1]-3*ksi[2]+2*ksi[0]*ksi[0]+4*ksi[0]*ksi[1]+4*ksi[0]*ksi[2]+2*ksi[1]*ksi[1]+4*ksi[1]*ksi[2]+2*ksi[2]*ksi[2];
  SF(1,0) = -ksi[0]+2*ksi[0]*ksi[0];
  SF(2,0) = -ksi[1]+2*ksi[1]*ksi[1];
  SF(3,0) = -ksi[2]+2*ksi[2]*ksi[2];
  SF(4,0) = 4*ksi[0]*(1-ksi[0]-ksi[1]-ksi[2]);
  SF(5,0) = 4*ksi[0]*ksi[1];
  SF(6,0) = 4*ksi[1]*(1-ksi[0]-ksi[1]-ksi[2]);
  SF(7,0) = 4*ksi[2]*(1-ksi[0]-ksi[1]-ksi[2]);
  SF(8,0) = 4*ksi[1]*ksi[2];
  SF(9,0) = 4*ksi[0]*ksi[2];
}

template<>
inline int QuadPoint<3>::mapDof(int elem, int dof, int iParent)
{
	throw 1;
}







template<>
void QuadPoint<2>::setSF()
{
  SF(0,0) = (1.0-ksi[0]-ksi[1])*(1.0-2.0*ksi[0]-2.0*ksi[1]);
  SF(1,0) = ksi[0]*(2.0*ksi[0]-1.0);
  SF(2,0) = ksi[1]*(2.0*ksi[1]-1.0);
  SF(3,0) = 4.0*ksi[0]*(1.0-ksi[0]-ksi[1]);
  SF(4,0) = 4.0*ksi[0]*ksi[1];
  SF(5,0) = 4.0*ksi[1]*(1.0-ksi[0]-ksi[1]);
}

template<>
void QuadPoint<2>::setGSF()
{
  GSF(0,0) = -3+4*ksi[0]+4*ksi[1];
  GSF(0,1) = 4*ksi[0]-1;
  GSF(0,2) = 0;
  GSF(0,3) = 4*(1-2*ksi[0]-ksi[1]);
  GSF(0,4) = 4*ksi[1];
  GSF(0,5) = -4*ksi[1];
  GSF(1,0) = -3+4*ksi[0]+4*ksi[1];
  GSF(1,1) = 0;
  GSF(1,2) = 4*ksi[1]-1;
  GSF(1,3) = -4*ksi[0];
  GSF(1,4) = 4*ksi[0];
  GSF(1,5) = 4*(1-ksi[0]-2*ksi[1]);
}

template<>
inline int QuadPoint<2>::mapDof(int elem, int dof, int iParent)
{
	return imesh->getVolumeNode(imesh->getFaceParentNo(elem,iParent),imesh->getFaceParentNodeNo(elem,iParent,dof));
}







template<>
void QuadPoint<1>::setSF()
{
  SF(0,0) = 0.5*ksi[0]*(ksi[0]-1);
  SF(1,0) = 0.5*ksi[0]*(ksi[0]+1);
  SF(2,0) = -ksi[0]*ksi[0]+1;
}

template<>
void QuadPoint<1>::setGSF()
{
  GSF(0,0) = ksi[0]-0.5;
  GSF(0,1) = ksi[0]+0.5;
  GSF(0,2) = -2*ksi[0];
}

template<>
int QuadPoint<1>::mapDof(int elem, int dof, int iParent)
{
	return imesh->getFaceNode(imesh->getEdgeParentNo(elem,iParent),imesh->getEdgeParentNodeNo(elem,iParent,dof));
}




template class QuadPoint<1>;
template class QuadPoint<2>;
template class QuadPoint<3>;


